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Abstract 

We study stochastic particle systems made up of heterogeneous units. We introduce a general 
framework suitable to analytically study this kind of systems and apply it to two particular models 
of interest in economy and epidemiology. We show that particle heterogeneity can enhance or 
decrease the size of the collective fluctuations depending on the system, and that it is possible to 
infer the degree and the form of the heterogeneity distribution in the system by measuring only 
global variables and their fluctuations. Our work shows that, in some cases, heterogeneity among 
the units composing a systems can be fully taken into account without losing analytical tractability 
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Mostjreal systems are made up of heterogeneous units. Whether considering a population 

m\ 



of cells 



2J , a group of people [3|, |4[ or an array of lasers |5|, |6[ (to name just a few examples) , 
one never finds two units which behave exactly in the same way. Despite this general 
fact, quantitative modeling most often assumes identical units, since this condition seems 
necessary for having analytically tractable models. Moreover, in the general framework of 
complexity science, systems very often can be modeled only at a stochastic level, since a 
complete knowledge of all the variables, the precise dynamics of the units and the interaction 
with the environment is not available. One way to include system heterogeneity is to consider 
that the interactions between the units are not homogeneous but mediated by some complex 
network, an approach that has attracted enormous attention in the last years 7|,|8j. An issue 
;hat has been less studied, beyond the role of particle heterogeneity in deterministic systems 
9r|l2|, is the heterogeneity in the behavior of the particles themselves in stochastic models. 
Some exceptions include the recent reference jl^ . where the authors analyze the effect of 
heterogeneous transition rates on consensus times in the voter model, and works considering 
the effect of a few "committed" individuals in this and related models [l4, [l^)]. In the 
context of statistical physics, the combined effects of stochasticity and heterogeneity have 
been considered, for example, in random-field Ising models and spin glasses 16|-18| or in 
diffusion in disordered media j3, Q] . We aim here at developing a general framework for the 
analytical study of stochastic systems made up of heterogeneous units, applicable beyond 
equilibrium models or Hamiltonian systems and suitable for a general class of complex 
systems of recent interest and at identifying some generic effects of particle heterogeneity 
on the macroscopic fluctuations. 

In this work we will show that the combined effect of stochasticity and heterogeneity 
can give rise to unexpected, non-trivial, results. While, based on naive arguments, one 
should conclude that global fluctuations increase in heterogeneous systems, we will show 
that in some systems of stochastic interacting particles fluctuations actually decrease with 
the degree of heterogeneity. Moreover, we will see that it is possible to infer the degree of 
particle heterogeneity (or "diversity") by measuring only global variables. This is an issue 
of great interest when one has access only to information at the macroscopic, population 
level, since it allows one to determine if heterogeneity is a relevant ingredient that needs to 
be included in the modeling. In this way, heterogeneity can be included when its presence is 
implied by the data and it does not enter as an extra free parameter. We will study first the 
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simple case of independent particles; then we will consider the general case of interacting 
particles and develop an approximated method of general validity to analytically study these 
systems; next, as a way of example, this method will be applied to two particular models of 
interest in economy and epidemiology. 

Our starting point is a stochastic description of a system composed by N non-identical 
units, which we call generically "particles" or "agents". Each particle is characterized by a 
constant parameter Aj (i — 1, . . . , N); the value of this parameter differs among the particles 
and it is the source of heterogeneity considered. Although there are more general ways of 



including heterogeneity, we will stick to this type of parametric heterogeneity [2l| because 
it is simple yet rather general. For simplicity, we assume that each particle can be in one of 
two possible states and define Si(t) = 0, 1 as the variable describing the state of particle i at 
time t (the two-states assumption will be relaxed later). The collective state of the system 
is given by the total number n(t) = ^2 i=1 s$(t) of particles in state 1. Sometimes, one does 
not have access to the individual dynamics and can only access experimentally the value 
of n{t). We are interested in the statistical properties of this global variable and how do 
they depend on the degree of heterogeneity in the system. We will often refer to n(t) as the 
macroscopic variable and to the Sj(t)'s as the microscopic ones. 



I. RESULTS 



A. Independent Particles 

We study first the case in which particles jump independently from state to 1 and 
vice-versa, schematically: 

rf r~ 

— 1, 1 — 0, (1) 

with rates that depend on the value of the heterogeneity parameter, rf = r ± (Aj). The 

probability Pi(t) for particle % to be in state 1 at time t obeys the linear rate equation 

-p- = —r^pi + rf(l — pi). In the case of constant rates, the solution is: Pi{t) = ^-(1 — 
dt 1 

e~ Tit ) + Pi(0)e~ Tit , with = rf + r~ . The results derived below apply equally if the rates 
depend on time or on the time that the particle has been in its current state (if the rate 
depends on the time a that the particle has been on its current state, the steady-state 



3 



probability of finding the particle at state 1 is pi >st = A +1 A - with = J °° dte fo dar t( a )) 

i ' i 

. Using particle independence and that the moments with respect to realizations of the 
stochastic process of the random variable Sj are given by (s^) = l k pi + fe (l — Pi) = Pi, one 
obtains that the average and variance of the global variable n are: 

N 

(*(*)> = £>(*) = (2) 



i=l 
N 



A<t)} = WO - rt f ) 2 ) = N (p^) ~ pW 2 ) > ( 3 ) 
i=i 

where the overline denotes an average over the population, g = 4 ^ 9i- ^ we consider a 
system where all particles are identical (i.e. have the same values for the internal parameter 
Aj = Aj,Vi,j), and keep the same average value (n(t)) for the global variable at time t, 
the variance would be of d [n(t)] = Np(t) (l - p(t)\ > a 2 [n(t)}. We obtain the somehow 
counterintuitive result that a system of heterogeneous independent particles displays smaller 
fluctuations in its collective variable than another system with identical particles. This effect 
is illustrated in figure [TJ The reduction in the variance of the collective variable is iV times 
the variance of Pi over the population: 

aUn(t)] - a 2 [n(t)} = N (pffi -^ff) , (4) 

which is of the same order, O(N), as the variance itself, giving a non- negligible correction. 

Reading the previous formula backwards, one realizes that the moments of the collective 
variable give information about the degree of heterogeneity in the system: 

-a (n(t))-(n(t)) 2 /N-a 2 [n(t)] 



p{tf -p{t) = ">' N LWi. (5) 

This expression is general, regardless the specific form in which is distributed over the 
population. Higher moments of the heterogeneity distribution are also related to higher 
moments of the collective variable. This allows to infer the skewness, kurtosis and higher 
order characteristics of the heterogeneity distribution by measuring only global variables 
and their fluctuations. In the Supplementary Information it is shown that an equivalent 
result is obtained generically for fc-state systems for k > 2. 

Besides the moments, one can derive the full probability distribution of the global vari- 
able. The generating function for the one particle probability distribution is g%(z) = 
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FIG. 1: Time series for the global variable n(t) of a system of identical (left panel) and hetero- 
geneous (right panel) particles, for a system of N = 100 particles. The parameters were set as 
rf = l,r~ = I/pi — 1, with pi = 1/2 in the case of identical particles (left panel) and pi chosen 
from a symmetric Beta distribution f(p) = [p(l — p)] Q_1 , with a = 0.05, being the sample 
mean and variance equal to p = 0.501, <J 2 \p] = 0.23, respectively. Note that the fluctuations of the 
average state are larger in the case of identical particles. 

^*. =0 z Sl P(sj) = 1 — Pi + PiZ and the generating function G(z) = YlH=QP{ n ) zn °f 
the sum of independent random variables is the product of their generating functions, 
G(z) = Ylf=i9i( z )- Expanding in powers of z we can obtain the probability distribution 

for n: P(n) = w!(A L n) , Eigg^ Ui=iPit Ue=n+i( l ~ PiJ» bein g 1 = •••,**) and S n the 
group of permutations of N elements. 

The model studied in this subsection, despite its simplicity, offers a reduced description 
of generic systems of non-interacting multi-stable units subject to fluctuations. The results 
obtained here are directly relevant if one is interested in the collective properties of one such 
system when the units are non-identical. One can reasonably argue that the independence 
property is too unrealistic for the study to be of any practical interest. We will be considering 
more complicated cases including non-independent units in the rest of the paper. However, 
this simple model presents in isolation a mechanism, spontaneous transitions, that can play 
a role in more complicated and relevant systems (we will see this later). The simplicity of 
the model allows us to understand the effect of heterogeneity in this mechanism, and will 
give us insight in the role of heterogeneity in the behavior of more complicated systems. 
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B. Two types of uncertainties 



We will now discuss the situation in which the particular values of the parameter of 
each particle are not known. This introduces an additional source of uncertainty. For 
simplicity we will focus on the 2-states independent particles system considered before, but 
the discussion applies as well to general systems of interacting particles. This discussion 
will also allow us to take a closer look at the results obtained and clarify their meaning and 
relevance in different settings. 

Often, one does not know the value of the parameter Aj of each individual particle, but 
has some idea about how this parameter is distributed on the population, perhaps its proba- 
bility distribution (obtained for example by measuring individual behavior in an equivalent 
system). Here, we will assume that the Aj's are independent and identically distributed 
random variables with a given probability density /(A). In this case, (n) and a 2 [n] are 
themselves random variables that, as shown above, depend on the particular realization of 
the Aj's. The expected values of these quantities are obtained by averaging Eqs. (l2f3j) over 
the distribution of the individual parameters: 

W0>=JVS*), <^W)\ = N (Sf) - pW) , (6) 

where the hat denotes an average with respect to /(A), g = f d\g(X)f(X). Again the 
variance is smaller than for a system of identical particles with the same mean value, namely, 

<j? d [n(t)} - ^Ht)\ =N(ffi- K*f) . 

If we average the probability P(n) over the distribution of parameters we obtain a simple 
form for the probability of the global variable n: 



P{n) = I dX x . . . dX N P(n\X u . . . , Atf)/(Ai) • • • f(X N ) 



n 



p n (l-pf~\ (7) 



a binomial distribution with parameter the average p = j dXp(X) f (X) over the distribution 
/. The variance of this distribution is 

a 2 [n(t)] tot = iv(p(tj-Sj 2 ), (8) 

equal to the variance one would obtain in a system of identical particles with the same 
average, Np(t), a result in apparent contradiction with (jfJJ). However, we should note that 



they refer to different things: Expression gives the average variance when the parameter 
values are given, so measuring the average uncertainty in n due to the stochastic nature of 
the process. (|S]) , in addition to the uncertainty coming from the stochasticity of the process, 
also includes the uncertainty on the parameter values. 

The two expressions are related by the law of total variance: 

(y 2 [n}tot = cr 2 \n] + a 2 [(n\X 1 , . . . , X N )] 



N(p-p 2 ) + a z 



E 



Pi 



N(p-f). (9) 



In cr 2 [(n|Ai, . . . , Ajy)], the variances are taken over the distribution of the A,'s. If we are 
considering a particular system, the temporal fluctuations (all the systems considered in 
this paper are ergodic, so we can think on averages over time or over the realization of the 
stochastic process interchangeably) in n will come only from the intrinsic stochasticity, and 
expressions f )3|6|) are the ones that measure it. Expressions flTfHj) are appropriate only if we 
are considering an ensemble of systems with a distribution of parameters and our different 
measurements may come from different systems in the ensemble. 



C. Formulation of the general method 

Let us now consider a general system of interacting heterogeneous particles. The stochas- 
tic description now starts from a master equation for the A^-particle probability distribution: 

dP(s 1 ,...,s N ) ^ r i 

j t = 2j.Ei - 1) [SiTi P(Si, . . . , S N )\ 

i=l 

N 

+ J2(E^ - 1) [(1 - ^k + P(si, . • • , s N )} , (10) 

with step operators defined as EfF(s%, Sj, ...Sjv) = F(s±, Sj + k, sn)- The transition 
rates rf might now depend on the state of any other particle (this is how interactions enter 
in the model). From Eq. lflO]) one can derive for the moments and correlations: 

d(si) 



dt 

d(siSj) 



(rt) ~ ((tT + rt)s t ) (11) 

(U -(qvsjad + irtsj + irpi) (12) 

with qij = + rj +rf + r~j and i ^ j in the second equation (recall that sf = si). In general, 
if the transition rates depend on the state variables s$, these equations are not closed since 



they involve higher order moments, and some approximation method is needed to proceed. 



Systematic expansions in 1/N, including van Kampen's f2-expansion 22], are not applicable, 
since variables Sj = 0, 1 are not extensive. In the following, we introduce an approximation 
suitable for the analytical treatment of systems of globally coupled heterogeneous particles. 

We assume that the m-particle correlations &ji,...,j m (t) = ($ji(t) ' ' 'fij m {t)) with 

Sj(t) = Sj(t) — {sj(t)) scale with system size as 

<rn,-,j m (t) = 0(iV- w/2 ), for j k ? j t . (13) 

Using this ansatz one can close the system of equations (I37|38l) for the mean values and the 
correlations. This is shown in the Supplementary Information for general transition rates of 
the form f( Sl /N,...,s N /N). 

While the resulting equations for the average values (sj(t)) coincide with the mean- field 
rate equations usually formulated in a phenomenological way[12l. 1 171]. our formulation allows 
us to compute the correlations and include, if needed, higher order corrections in a systematic 
way. 

Assumption (139|) can be justified noting that it is consistent with ^ ^ j m a h,...,j m = 
((n-(n)) m ) = 0(N m / 2 ) which follows from van Kampen's splitting of the global variable n = 
N<p + iV 1 / 2 ^, with (ft deterministic and £ stochastic. Details are given in the Supplementary 
Information. The global variable n is extensive and it is expected to follow van Kampen's 
ansatz in many cases of interest. Note, however, that since there is not a closed description 
for the macroscopic variable n, one can not use van Kampen's expansion, and our approach 
extends the implications of this splitting of the macroscopic variable to the correlations of the 
microscopic state variables. For simplicity, we have focused on 2-states systems and assumed 
a constant number of particles. Systems with k states are also expected to follow ansatz (13TJ1) , 
since the scaling of the global variable is not limited to 2-sates systems. The case of variable, 
but bounded, number of particles can be included straightforwardly by considering an extra 
state. The unbounded case can also be considered performing an appropriate limit. If the 
system has some spatial structure, the ansatz (139|) is not expected to be valid, and some 
decay of the correlations with the distance is expected instead; this interesting situation is 
left for future work. 

We will proceed by applying the presented method to analyze the role of heterogeneity 
in two models previously considered in the literature that apply to contexts in which the 
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assumption of identical agents can hardly be justified: stock markets and disease spreading. 
We will focus on the steady-state properties of both models, skipping transient dynamics. 



D. Application to Kirman Model 



Kirman's model |23j was proposed to study herding behavior in the context of stock 
markets and collective dynamics on ant colonies. In the stock market context, agent i 
can be in two possible states (e.g. = "pessimistic" -with regard to future market price- 
and 1 = "optimistic" ) and it can switch from one to the other through two mechanisms: 
spontaneous transitions at a rate e, and induced transitions at a rate A -1 £\ ^jO- ~~ ^s t , Sj ), 
being A,- the "influence" of agent j on other agents. The case e = corresponds to the voter 



model 



24[ | . In the original formulation, all agents have the same influence, i.e. A, = A-,-, Vi, ?'. 



We generalize the model allowing the parameter Aj to vary between agents. In 25], the 
effect of heterogeneity was explored numerically, but not in a systematic way. 

This model is interesting for us because it incorporates in a simple way two basic processes: 
spontaneous transitions and induced transitions. As we will see, due to its simplicity, a full 
analytical treatment is possible that will, in turn, allow us to obtain a deeper insight into 
the general effect of heterogeneity in systems of interacting particles. 

The master equation for the process is of the form ( JTUI) . with rates given by: 

r+ = e + N- 1 A***, rf = e + A^ 1 £ A fc (l - s k ) (14) 
k k 

From (1371) the averages and correlations obey: 

= e _ (2e + A)( Sl )+A- 1 ^A fe ( Sfe ), (15) 

k 

da id 



2(2e + \)a id + A" 1 ^ A fc (a ijk + a jjk ) (16) 

k 

for i ^ j and a iti = (sj)(l — (sj)). Note that, due to the particular form of the rates, 

these equations do not involve higher-order moments. This is a simplifying feature of this 

model that allows one to obtain exact expressions. The first equation leads to a steady state 

value (n) st — y (a property that comes from the symmetry <H- 1). Using the relation 

a2 [ n ] = Y2i j a i,j we obtain (see Supplementary Information) that the variance in the steady 

state is: _ 

M T 2A(l-iV~ 1 ) (A - 3 + 2A- 1 U" 

1 + — ^ + ~ ; = — (17) 



2r 1 N 



4e + A 2e + A 
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FIG. 2: Variance of the number of agents in state 1 as a function of the variance of the influence 
parameter A in Kirman's model with distributed influence. Numerical simulations (symbols) and 
theoretical analysis (lines), Eq. (|53p . for different number of agents N and e = 0.01. Aj are indepen- 
dent random variables distributed according to a lognormal or a gamma distribution with mean 
A = 0.5 and variance cr?. The results have been averaged over 2 x 10 4 for N = 50 and 10 4 for 
N = 100 realizations of the distribution of parameters. 



N 



with Ai = N{4e+ j )+2X% ■ The leading-order term, a 2 [n] , 



1 + A + ^L 



2e 2e(4e+A) 



+ O(N ), With 



a 2 [A] = A 2 — A 2 , can also be readily obtained using the ansatz (|39|) . Note that the presence 
of heterogeneity increases collective fluctuations. In Figj2] we compare expression f l53|) with 
results coming from numerical simulations. 

In this case, the knowledge of (n) st and ofjra] alone does not allow to infer the degree 
of heterogeneity present in the system, unless one knows from other sources A and e and 
it is not possible to conclude whether the observed fluctuations have a contribution due to 
the heterogeneity of the agents. However, the steady-state correlation function K[n](t) = 
(n(t)n(0)) st — (n)l t , does include a term that allows to infer the possible heterogeneity. 
K[n](t) is obtained integrating Eq. (j45il and performing the appropriate conditional averages 
(see Supplementary Information): 

K[n](t) = (4M - u) e-^ +J » + ue- 2 '\ (18) 

with u = ^(f^/jv) ( a stN — N/4). The departure from a pure exponential decay signals the 
presence of heterogeneity (for identical particles u = a^ t [n]). Fitting this expression to data 
one can obtain ofjn] and the parameters e, A. Then, the use of expression f )53|) would yield 
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FIG. 3: Correlation function (in log-linear scale) for Kirman's model with distributed influence. 
Results coming from numerical simulations (symbols) and theory (Eq. (|18p . solid lines). Note that 
when heterogeneity is present (er^ = 1.5) the correlation function departs from purely exponential 
decay (displayed as a dashed line). Data for o\ = have been moved up 5.5 units vertically for 
better visualization. Parameters values are e = 0.01, N = 100. Aj are independent random variable 
distributed according to a gamma with mean A = 0.5 and variance, a\, indicated in the figure. A 
simple fit of expression (fTBj) to the a\ = 1.5 data gives A = 0.50, e = 0.0099. 

a 2 [A] . In Figj3] we show that the numerical simulations indeed support the existence of two 
exponential decays for the correlation function. 

1. Other ways to introduce heterogeneity 

Interestingly, other ways to introduce heterogeneity in the system have different effects: 
-If the heterogeneity is introduced in the spontaneous transition rate, e — > Cj, making some 
particles more prone to spontaneous transitions that others (but keeping Aj = A, Vj to 
isolate effects), collective fluctuations again increase with respect to the case of identical 
particles. 

-Next, we can assume that the rate of induced change is different for different agents, even 
if all have the same influence. Measuring this difference in "susceptibility" (to induced 
change) with a parameter Ui, we would have that the rate of induced change in agent i is 
uii ^ • Xj(l — d s . ;Sj )/N. The effect of heterogeneity in uj; t (keeping again Aj = A, Vj) is that 
the collective fluctuations decrease with the degree of heterogeneity in the susceptibility Wj. 
-Setting some heterogeneous preference for the states among the particles, i.e. making ef, 
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the spontaneous rate from to 1 of particle i, different from e~, the spontaneous rate from 
1 to of the same particle, decreases global fluctuations. In order to vary the preference 
for one state keeping constant the global "intrinsic noise" of this particle (note that the 
correlation time of particle i, when isolated, is given by e+ + e 4 ~), we set = 2e — e~ and 
generate as i.i.d. random variables with a distribution with support contained in the 
interval [0, 2e] . Exact explicit expressions for the first moments of the global variables are 
(see Supplementary Information): 



In figure (J4j) the exact expressions ( JT9| 154)) are compared with results coming from numerical 
simulations. 

In this case, the correlation function decays exponentially, 



independently of the degree of heterogeneity, so this form of heterogeneity cannot be inferred 
by measuring the correlation function. Numerical simulations confirm this result. 

2. Intuitive explanation of the effects of heterogeneity 

We have seen that heterogeneity can have an ambivalent effect over the size of the fluc- 
tuations, depending on the particular form it appears. We now provide intuitive arguments 
to understand these different effects. 

When the influence parameter, Aj, varies from one unit to the other, there will be some 
largely influential agents and others with little influence. In the limit of very large hetero- 
geneity we can think of a situation with a single agent with an extremely large influence 
and the others having a negligible one (we are keeping a constant average influence). In this 
case, the highly influential agent drifts from one state to the other, essentially independently 
(since other agents have negligible influence), but, due to its large influence, all the agents 
are attracted to its current state. In this "follow the leader" regime, we obtain macroscopic 
transitions from one state to the other, corresponding to very large global fluctuations. 




(19) 



(20) 



K[n]{t) = a 2 [n] st e 



(21) 
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FIG. 4: Variance and average of the number of agents in state 1 as a function of the variance of 
the spontaneous transition rate to state 1, e + , in Kirman's model. Results coming from numerical 
simulations (symbols) and theoretical analysis (solid lines, Eqs. (jl9l I54p ). for N = 50 agents, Aj = 
A = 0.5 and e+ + e_ = 2e = 0.4. ef are independent random variables distributed according to a 
symmetric beta distribution in the interval (0, e) with mean e + = 0.2 and variance, a 2 + . 

The situation is the opposite for a non-identical susceptibility parameter Ui where global 
fluctuations decrease as the diversity is increased. Again, we can understand this in the 
limit of very large heterogeneity where a single agent (or a small number of them) has large 
susceptibility while all the others have a negligible one (in order to keep average susceptibility 
constant). Then, agents with small susceptibility change essentially independently, in an 
uncorrelated fashion, resulting in low global fluctuations (note that in order to have large 
global fluctuations, the fluctuations in the state of the single agents should be correlated). 

In the case of diverse spontaneous transition rates, e^, global fluctuations increase with 
the degree of heterogeneity. In the limit of large heterogeneity, we would have a small 
number of agents with very large spontaneous transition rate, whose state would fluctuate 
in an uncorrelated fashion, and a large number of agents with low spontaneous transition 
rate, that essentially would only change state through induced transitions, giving rise to 
correlated fluctuations, resulting in large variance for the global variable. 

In the case in which agents display an intrinsic heterogeneous preference for one of the 
two states, the global fluctuations decrease with heterogeneity degree. We saw this already 
in the first section for non-interacting agents. Here we see the same effect, suggesting that 
the phenomenon is robust and still plays a role when interaction is added. 

The coexistence of a small number of agents with large value for the parameter and a 
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large number with a small value assumed in the previous arguments arises from the fact 
that an unbounded (from above) distribution of positive-defined parameters (e.g. rates) is 
skewed. However, all the effects of diversity commented are still present if the distribution is 
symmetric. In this case, nevertheless, the maximum degree of heterogeneity (for a constant 
mean value) is bounded, sometimes greatly limiting the maximum possible value of diversity. 
For symmetric distributions, a simple explanation is not so clear, but an asymmetry in the 
effect of increasing and decreasing the value of the parameter seems to be at the heart of 
the phenomenon. 



E. Application to the SIS disease spreading model 

The previous example could be treated exactly because, due to symmetry, the interaction, 
non-linear terms, cancel out in the equations for the moments. In general, however, this is not 
the case, and the analytical treatment is more involved. Here we consider an example of such 
case. The stochastic susceptible- infected-susceptible (SIS) model is a paradigmatic model for 
the study of spreading of infectious disease {26] as well as the diffusion of innovation |11| and 
other types of social influence. Despite its simplicity, it captures interesting phenomenology. 
The process is schematically described by: 

\ 3 /N 7 

s(i) + -> m + m, iu) -> su), s(j) -> I®, (22) 

where S(i) (resp. denotes agent i being susceptible (resp. infected). There are 3 

basic elementary processes: (i) infected agent j infects susceptible agent i at a rate Xj/N, 
being Xj the infectivity parameter of agent j; (ii) infected agent j becomes susceptible a rate 
7; (iii) susceptible agent j gets infected spontaneously (due to interactions with agents not 
considered in the system or other causes) at a rate e. This corresponds to the SIS model with 
spontaneous contagions and distributed infectivity. In the absence of spontaneous infections 
e = 0, the system has a trivial steady state with zero infected agents. With e ^ the 
system has a non-trivial steady state whose properties we analyze in the following. As in 
the previous case, heterogeneity could appear in any parameter of the agents (for example, 
in the recovery rate, in a "susceptibility" parameter, etc.). 

We study first the case in which only the infectivity, Aj, can vary from agent to agent. The 



effect of heterogeneity in the deterministic version of related models was studied recently 12]. 
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FIG. 5: Average and variance of the number of infected agents in the SIS model as a function of 
the variance of the infectivity. Numerical simulations (symbols) and theoretical prediction to first 
order (lines). Parameters values are e = 0.01, 7 = 1, N = 200, A = 0.5. 

The master equation is of the form (jTOl with rates rf = e + ^jj^-,r^ = 7. Equations (1371 - 
[38]) for the first moments can be closed in the steady state, using our main ansatz, to obtain 
explicit formulas for (n) st and cr 2 [n] st to any desired order in A r_1 . In this case, however, 
the expressions are rather cumbersome and we skip them here. The results are plotted 
in figure ([5]), where we compare the approximation to order 0(A r_1 ) with results coming 
from numerical simulations, showing good agreement. Here both the average value and the 
variance are modified by the presence of heterogeneity (the dependence of the average is, 
however, only in second order in 1/N, almost unnoticeable in the figure). As in the Kirman 
model, the size of the fluctuations increase markedly with the amount of heterogeneity in 
the "influence" (now influence to infecting others) of the agents. 

In this case, other ways to introduce heterogeneity also have different effects. When 
heterogeneity appears in the recovery rate 7, the mean number of infected agent increases, 
with a moderate effect over the variance (resulting in smaller relative fluctuations). Het- 
erogeneity in the susceptibility to infection (which would be introduced with the change 
rf = e + ^ ^j^- — > e + Ui ^j^-, with U{ distributed over the population) decreases 
the fluctuations, with little effect over the mean value. Heterogeneity in the spontaneous 
infection rate e has almost no effect. In a real situation, one expects to find heterogeneity si- 
multaneously in several of the parameters defining the model. When heterogeneity is present 
both in the infectivity and in the susceptibility, the effects of both types of heterogeneity 
essentially add up, with the size of the fluctuations increasing with the heterogeneity in the 
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infectivity for a given level of heterogeneity in the susceptibility and fluctuations decreasing 
with the level of heterogeneity in the susceptibility for a given level of heterogeneity in the 
infectivity. The effects of heterogeneity in the infectivity and in the susceptibility are equiv- 
alent to those found in the Kirman model, and can be intuitively understood in the same 
terms. Heterogeneity in the recovery rate is similar to assigning an heterogeneous preference 
for the state (recovery) and its effect in the (relative) fluctuations is again the same as 
that in the case of the Kirman model. This suggests that the effects of the heterogeneity 
found are generic and can be useful to understand the behavior of other systems. 

II. DISCUSSION 

In this work, we have analyzed the combined effect of stochasticity and heterogeneity in 
interacting-particle systems. We have presented a formulation of the problem in terms of 
master equations for the individual units, but extracted conclusions about the fluctuations 
of collective variables. We have developed an approximation suitable for the analytical 
study of this general type of systems. We have shown that the heterogeneity can have an 
ambivalent effect on the fluctuations, enhancing or decreasing them depending on the form 
of the system and the way heterogeneity is introduced. In the case of independent particles, 
heterogeneity in the parameters always decreases the size of the global fluctuations. We 
have also demonstrated that it is possible to obtain precise information about the degree 
and the form of the heterogeneity present in the system by measuring only global variables 
and their fluctuations, provided that the underlying dynamical equations are known. In 
this way stochastic modeling allows to obtain information not accessible from a purely 
deterministic approach. We have also demonstrated that, in some cases, one can account 
for the heterogeneity of the particles without losing analytical tractability. 

Heterogeneity among the constituent units of a system is a very generic feature, present 
in many different contexts and this work provides a framework for the systematic study 
of the effect of heterogeneity in stochastic systems, having thus a wide range of potential 
applicability. More research in this direction would be welcomed. 
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III. METHODS 



We have developed and used analytical tools based on an extension of van Kampen's 
ansatz on the relative weight of the fluctuations compared to the mean value, suitable for 
systems with particle heterogeneity. We have included the details of this method in the 
Supplementary Information. In some cases, and in order to compare with the analytical 
expressions, we have generated data from numerical simulations using a particular form of 
Gillespie's algorithm |27j that takes into account the heterogeneity in the population. We 
now explain this algorithm using, for the sake of concreteness, the specific case of the Kirman 
model with distributed susceptibility. 

The parameters of the system are: the spontaneous transition rate e, the influence pa- 
rameter A, the susceptibility parameter of each agent Ui, and the total number of agents N. 
In this case, the influence parameter can be reabsorbed rescaling u i: so we set A = 1 without 
loss of generality. The variables of the system are the state of each agent s» = 0, 1. We will 
also use the total number of agents in state 1, n = ^2 Sj, the total susceptibility of agents in 
state 1, ft — ^2 SiUi, and the average susceptibility to = 

At any given instant, two events can happen: 
(i) An agent in state 1 changes to state 0. This can happen due to a spontaneous transition, 
at a total rate ne, or due to an induced transition, at a total rate £1 . 
(i) An agent in state changes to state 1. This can happen due to a spontaneous transition, 
at a total rate (N — n)e, or due to an induced transition, at a total rate (NoJ — 

According to the Gillespie method, that considers the continuous-time process, the time at 
which the next transition will take place is exponentially distributed, with average the inverse 
of the total rate. The probability that a given transition is realized is proportional to its rate. 
If the realized transition is a spontaneous one, the agent that actually undergoes it is selected 
at random (since, in this case, they all have the same rate e). If the transition is induced, the 
agent that undergoes it is selected with probability proportional to its susceptibility. It can 
be easily seen that this principles lead to an exact (up to numerical precision) simulation of 
sample paths of the stochastic process 

The algorithm, then, proceeds as follows: 
(0) Evaluate the total number of particles in state 1, n, and the total susceptibility of 
particles in state 1, Q. 
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(1) Evaluate the total transition rate r = eN + fi^p + (Nco - . 

(2) Generate the time for the next reaction, t n , as an exponential random variable with 
average 1/r. This can be done by setting t n = — Mn£7, with U a uniform random variable 
in the range (0, 1). 

(3) Select which reaction takes place. For this, generate a uniform random variable, g, in 
the range (0, r). 

-If g < ne, the transition will be a spontaneous transition from 1 to 0; Select an agent, j, 
at random among those at state 1. Set n — n — 1, Q = Q — u>j. 

-If ne < g < Ne, the transition will be a spontaneous transition form to 1; Select an 
agent, j, at random among those at state 0. Set n = n + 1, Q = Q + Uj. 

-If iVe < g < Ne + n ( j ^~ Tt ) ; the transition will be an induced transition from 1 to 0; 
Select an agent, j, among those at state 1 with probability proportional to the value of its 
susceptibility parameter Uj. Set n = n — 1, Q — Q — Uj. 

-If g > Ne + Qiitpl) ; the transition an induced transition from to 1; Select an agent, 
j, among those at state with probability proportional to the value of its susceptibility 
parameter Uj. Set n — n + 1, Q = Q + Uj. 

(4) Set t = t + t n . Go to (1). 
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Supplementary Information 



A. M-states system 

We consider here the case in which each particle can be in one of M (instead of 2) possible 
states. We will show that the results obtained in the main text for 2— state systems also 
hold in this more general case. 

We label the states with the subscript a — 0, 1, . . . , M — 1, so in this case the variable 
describing the state of particle % can take M possible values, = 0, . . . , M — 1 (we start the 
labeling from to be consistent with the previous case, that would correspond to M — 2). 
Let pi(Xi,a, t) the probability that particle i, with heterogeneity parameter Aj, be on state 
a. It satisfies the evolution equation: 

^ (A i; M) = E^AOpA./M), (23) 



dt fi 



with A a> p a general transition matrix (satisfying J2^=q A lj0t = 0), that may depend in 
principle on time and on the time that the particle has been on its current state. To isolate 
the role of parameter heterogeneity, we assume that the initial condition is the same for all 
the particles (or that the initial condition is determined by the value of Aj) such that the 
solution pi(Xi,a, t) = p(Xi,a,t) is the same for all particles sharing the same value of the 
parameter. The macroscopic state of the system will be described by the set of variables 
n a = J2f =1 8 a , Sl , that is, the number of particles in each state. The averages and variances 
of this variables are given by: 

N 

(n a (t)) = ^p(A;,cM) (24) 



N 



An*(t)} = J2 [p(Ka,t)-p(\ t ,a,t) 2 ] . (25) 

i=l 

This variance is again smaller that tat of a system of identical particles with same average, 
the difference given by: 



v 2 [n*(t)] id - a 2 [n a (t)} = Np(a, tf - p(a, t) , (26) 

a result exactly analogous to the one obtained in the previous case. The heterogeneity 
among the particles on the probability of occupation of level a can be derived from the first 
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moments of the occupation number of the level: 

-2 (n a ) - (n a ) 2 /N - cr 2 [n a } 



p(a,t) 2 - p(a,t) 



N 



(27) 



Note that, when focusing on the number of particles on state a, the system effectively 
reduces to a 2— level one, with states a and no- a, so the results of the previous section can 
be translated directly. 

A different and some times relevant question can be considered when the labeling of the 
states is such that the order is well defined (for example each state corresponds to an energy 
level or a distance from a reference). Then the average state is meaningful and we can study 
its statistical properties. Below we show that the variance of this mean level is again always 
smaller if heterogeneity is present. 

The average state of the system is given by L = ^2^=0 a It is a random variable 
whose average and variance are given by: 



(L) 



M-l 

E 



a- 



o=0 

M-l 



N 



M-l N 

EE 

o=0 i=l 



p(Xi,a) 



a- 



N 



N 



a,P=Q i=l 

We have used p(\i, a) = (S a ,si) and (n a np) 



(28) 

'M-l M-l 

o,/3=0 



o=0 



J2i,j=l(^a,SiSp, Sj ) 



{n a )(np) + 



Ylii=i[^a,PPi a i ^i) ~ P( a i ^i)p(Pi \)]- A system of identical particles that had the same aver- 



age occupation of the different levels i.e. pid(\i,a) = jf~^2f=iP(^ji a i) = Vi, a, would 
have and average and variance of the mean level given by: 

M-l 



Q=0 

M-l 



N 



a 2 [L] 



M-l 



id 



(n a ) (np) 
N N 



(30) 
(31) 



o=0 o,/3=0 

We now define g(\i) = ocp(Xi, a) (the average level of particle i), and note that the 
first terms in the right-hand side of ( )29l) and (13T1) are equal, while the second terms can be 
written as: 

N M-l AT 

(32) 



^ N M-l N 

mYY ap(Ka)MKP) = 9{\? ! = T7F, 



7V2 ^ /-^ -rw-irrwr, ^ 2 

i=l o,/3=0 i=l 



M-l 



a,8=0 



(wq) {np) J_ 
N N ~ N 



1 N 

nY^ 

i=l 



N" 

2 



N 



9 2 , 



(33) 
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which implies that a 2 [L]id > er 2 [X], i.e. the variance of the mean level is always smaller in 
a system of heterogeneous particles, the difference with respect to the case of identical ones 
being: 



M-l 



a,0=O 



,1=1 t,j=l 



> 0. 



(34) 

The correction to the variance in this case scales as 1 /N, but again is of the same order as the 
variance itself, indicating a non-negligible correction. In this case to derive the heterogeneity 
of g(Xi) over the population one needs to know the average occupation level of each state 
(n a ) and use: 

¥-f = Yl °?M/N - (L) 2 - Na 2 [L\. (35) 

a 

This can be written in terms of the variance of L in an equivalent system of identical particles, 
<7 2 [L] id . If this is known, one can directly use 

7-f = N(a 2 [L] id -a 2 [L]). (36) 

Note that, contrary to the two- level case, now the value of (L) does not determine <r 2 [L] id . 



B. Intuitive origin of the decrease of fluctuations for independent units 

We have shown that a system of independent heterogeneous particles has smaller fluctu- 
ations for the collective variable than an equivalent system of identical ones. The origin of 
this result is the following (for simplicity we refer to the case of 2-state system): 

The average of the global variable is determined by the concentration of the states of the 
particles around state 1 ((n) = J2i( s 'd)- The fluctuations (measured by the variance) of 
the global variable are determined by the stochastic fluctuations of the individual particles 
alone (a 2 [n] = X^°" 2 [ S «L since the particles are independent). 

In a system of heterogeneous particles, the dispersion of the states of the particles is due to 
the heterogeneity (some prefer to be around sate 0, others prefer to be around sate 1) plus 
their intrinsic stochasticity. In a system of identical particles, the dispersion comes from 
the stochasticity alone, so for a system of identical particles to have the same concentration 
in the states of the particles (global average) than a heterogeneous system, the intrinsic 
stochasticity has to be larger. This will give rise to larger fluctuations for the global variable. 
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In particular, any given rational value of ^ = ^ can be obtained with zero fluctuations, 
taking A particles that are always at state 1 and B — A particles that are always at state 0. 

This explanation is illustrated in figure (EJ). In the identical-particles system both particles 
fluctuate between 1 and 0. In the heterogeneous case, one particle spends most of the time 
at 1 and the other spends most of the time at 0. The probability of finding a given particle 
at 1 is the same in both cases (1/2) but in the heterogeneous case most of the time there 
is one particle at 1 and one particle at 0, resulting on a value of the average state most 
often equal to 1/2, and so with smaller fluctuations. The situation is similar for a larger 




0.8 

0.5 - 



FIG. 6: Time series of a system of two identical (upper panel) and heterogeneous (lower panel) 
particles, together with the corresponding average state. Note that the fluctuations of the average 
state are more pronounced in the case of the identical particles. 

number of particles, as shown in figure 1 of the main text. An analogous picture emerges 
when one considers more that 2 states. Note that in every case we compare a system 
of heterogeneous particles with another of identical ones that has the same one-particle 
distribution i.e. Pi(a) ic [ = V. p ^ a ^ , Vi,a. 



C. Justification of the Ansatz 



The general evolution equations for the first moments are of the form: 

d(si) 



dt 

djsjSj) 

dt 



(rt)-((r7 + r?)*<), (37) 
-(qijS^ + irts^ + irpi). (38) 
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Our main ansatz is that the m-particle correlations o'j 1 ,...,j m {t) = (Sj^t) ■ ■ • Sj m (t)) with 
8j(t) = Sj(t) — (rij{t)) scale with system size as 

°n,., jm (t) = 0(N- m / 2 ), farj^jt. (39) 

We first show how the ansatz (1391) allows to close the system (|3"T| 1381) . 
We assume that functional dependence of the rates on the sate variables is of the form 
f(si/N,...,sjsr/N). This includes, for example, rates of the form f(%2 A^s^/iV) like the 
ones used in the examples analyzed. We further assume that the rates can be expanded as 
a power series: 

N ^ N I N 

ii=l ii,i2=l ii,...,ifc=l 

(40) 

JV 

There are N k terms in the fc'th summand, , giving a total contribution of order 

h,...,i k =l 

0(N°). The terms in the right hand side of (1371) are of the form: 

(fii • • • Q _ (gn + (Q) ; ; ; (K + (Q)) _ ^ g<g)^ _ ^ Q(iV^ /2 ) 

fc! jfe! Z^i!(jfe_i)! /!(£;_/)!' 1 J 

where 5 l corresponds to a term of the form (5j 1 (t) ■ ■ ■ 5j l (t)) , (s) k ~ l corresponds to 
(s^) • • • (si k _ t ) and the last equality holds due to our ansatz. We see that the dominant 
terms are those with 1 = 0, which correspond to products of mean values of the form 
(s^) ■ • • (si k ). We conclude that the ansatz allows to do the substitution (s^ . . . s ik ) — > 
(sii) " " " ( s i k ) + 0(N~ 1 / 2 ) in the evolution equations for the mean values. 
The evolution equations for the correlations read: 
da* 



• -j 



((rr + rf) Sl 5 3 ) + ((r7 + r+) 8j St) + (r+Sj) + (r+Sj). (42) 



dt 

In this case, the terms are of the form (s^ . . . Si k 5 r ) = ((8^ + (s^)) . . . (Si h + (si h ))5 r ) with 
r = Due to the presence of 5 S , the term in which only averages appears vanishes. 

Reasoning as before, we see that the dominant terms are those proportional to (Ti hS , while 
those proportional to higher-order correlations can be neglected. In this case, the ansatz 

k 

allows to do the substitution (s^ . . . s ik 5 r ) — > (s^) ■ ■ ■ (s ik ) v + 0(iV _3//2 ). In this 

1=1 ^ 

way, the evolution equation for the correlations depend, at first order, only on averages and 
correlations and not on higher order moments. 
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The validity of the ansatz (139|) itself can be established a posteriori by checking that 
the results obtained using the ansatz are consistent with it. In this section, we will link its 



validity with the well-known van Kampen's ansatz 22| that is the basis for the systematic 
system-size expansion. 

Van Kampen's ansatz consists on assuming that the variable of interest has a deterministic 
part of order Q plus a stochastic part of order fi 1 / 2 , i.e. n = Vtcj){t) + fi 1 / 2 ^, where Vt is a 
parameter of the system that controls the relative size of the changes due to elementary 
processes, typically the system size. 

In our system the role of the parameter Q is played by the total number of particles N. 
As briefly stated in the main text, we cannot expect that the single-particle variables that 
we are considering obey van Kampen's ansatz, since they are not extensive. Our variables 
Si — 0,1 have a deterministic and stochastic part that are both of order zero respect to iV 
(note that cr 2 [sj] = — (s»))). However, the macroscopic variable n = J^Sj is indeed 

extensive and we can expect that it will follow van Kampen's ansatz: n = N<fi(t) + N 1 ^ 2 ^. 
This implies that the m-th central moment of n will scale as 7V m / 2 , i.e: 

<(n-<n» m >= £ °r j 0{N m ' 2 ). (43) 

Now, assuming that crj l! ...j m = fm{N)d-j lt _ t j m for j k ^ ji, with aj lt _ t j m independent of N 
i.e. the m-particle correlations are all or the same order in N, so that X^i^w +j m &ji,-,jm 
scales as N m (note that there are of the order of N m terms in the sum), we obtain our main 
ansatz, <7 3 \,...,j m = 0(N~ m ^ 2 ) for jk ^ ji- We have only considered terms with jk ^ ji in 
the sum (H3l) ; terms with repeated sub- indexes can be expressed as lower order ones. For 
example, if the index j\ is present k times, and the others are all different, we find: 



a jl,jl,---jl,j2,---j m -k+l (i S jl ( S jl)) ^32 ■ ■ ■ fijk-k+l) 

= ^ 2 ,..j m - fc+ i(-(sj 1 )) fc + (Sj, ■ ■ ■ s jm _ k+1 ( . (-(sj^Ysj, 

i=0 w 



) (44) 



32,— Jm 



as can be see expanding (sj 1 — (sj 1 )) k and keeping in mind that s 2 = Sj. The number of such 
terms in the sum ( H3"|) is 0(N m ~ k+1 ), so they give smaller contribution that terms with all 
sub- indexes different. Proceeding order by order from k = 1, we see that our main ansatz 
(EHD follows from (03}. 
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We point out that in systems of heterogeneous particles we do not have a closed description 
for the global, extensive, variable n so van Kampen's expansion cannot be used. Instead we 
derive the implications of van Kampen's ansatz over the correlations of the microscopic vari- 
ables, is a simple and convenient expression that in general allows to close the equation 
for the moments (I37|I38I) . Often, however it is not necessary, and a weaker condition of the 
form (14*31) . that directly follows from van Kampen's ansatz without further assumptions, is 
sufficient. 

Van Kampen's ansatz is generally valid when the macroscopic equations have a single 
attracting fixed point, when the system displays small fluctuations around the macroscopic 
state. The general method explained here is expected to be valid under similar conditions. 
An interesting topic for future research will be whether a system that has a single attracting 
fixed point in the absence of diversity always maintains this globally stable state when 
diversity is present, and whether a system that does not posses this globally stable fixed 
point can acquire it when diversity is added. 



D. Details of the calculation for the Kirman model 



In the Kirman model with distributed influence, the averages and correlations obey: 
d{Sl) ' = e-(2e + A)( Sj ) + iV- 1 ^A fe ( Sfe ), 

k 

" - -2(2e + X)a id + N' 1 ^ A fc {a i%k + a 



dt 
dai 



(45) 



dt 



e + a+{X-2a)(si) -2^ 



N 



(46) 



with a = Afc j^ . Note that, due to the particular form of the rates, these equations are 
indeed closed. The first equation leads to a steady state value (si) st = |, which implies 



(n) 



St 



Y (a property that comes from the symmetry -H- 1). (j46p is a linear system 



of equations for the correlations. The steady state correlations can always be obtained by 
inverting the matrix that gives the couplings. Obtaining a closed expression for a 2 [n] in 
terms of the moments of A is, however, not straightforward. From ( |4*6|) . we see that in the 
steady state: 



'■j 



E k A fc ^4^ + &u [e + A/2 - 2 Z k 
2(2e + A) 



(47) 
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from where we obtain 



a 2 [n] = £ 



a 



N(e + A/2) + 2C(1 - 1/JV) 



2(2e + A) 

with C = £\ . XjUij- Multiplying ( T4"T|) by Aj and summing over j, one obtains: 

d(l-2/N) + (e + X/2)XN 



(48) 



C 



4e + A 



(49) 



where d = ^ i ■ AjA^er^-. This is obtained again multiplying (T4T|) by AjAj and summing over 



d 



(e + A/2)(A 2 )iV-2e/iV 
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(50) 



where e = Ylij ^i^j a i,j- Using the ansatz a^j = 0(N x ) we see that the last term of ( 1501) i 
0(N°) (while the other are of O(N)), so to the first order we obtain : 



N 

T 



l 



A 



a 2 [A] 



2e 2e (4e + A) 



+ 0(JV°), 



(51) 



with a 2 [A] = A 2 - A . We have seen how the application of the ansatz ( 1391) allows one to 
obtain closed expression for the global average and variance. Interestingly, in this particular 
example, it is possible to include all higher order terms to obtain an exact expression for d 
(which gives the exact expression for a 2 [n] trough (I49|48p ). details are given in the appendix: 

fc T2Tfc jV(e + A/2)^g^ 

(52) 



d 



N(e + A/2) J21 



k=0 \ N(4e+\) 



46 + A-Er=o 



_ JV(4e+A) 

The second equality holds as long as lim 



4e + 



2A 2 



iV(4e+A)+2A 
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N{Ae+\) 



(note that 
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J2a k X k + 2 = A 2 ^a fc A* 



fe+l \k+l 



X 2 



k=0 

(A+4e)jV 



k=0 



, 1 - a k+1 X 
1 — aX 



JV(4e + A) 

A sufficient condition for this is A, < 



N. When the Aj's are i.i.d. random variables, the probability that 



this condition is satisfied approaches one as N grows. This condition is actually necessary 
and sufficient for the first equality in ( |52l) to hold (see appendix). 
We finally obtain the following exact expression for the variance: 



2r 1 N 



A 2 



l+ 2A(1 ~!_ /iV) +(iV-3 + 2/iV)- V(1 " + " )+2A 



4e + A 



2e + 



A 2 



jV(4e+A)+2A 



(53) 
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We see from (|53|) that higher order corrections to a 2 [n] depend on higher order moments 
of the distribution of A over the population. 

Expressions (IBTj refer to the variance of n in a population with given values for the pa- 
rameters of each agent, A«, so the averages are population averages i.e. g(X) = Yli=i 9(\)/N- 
In the case that the parameters of the agents are random variables, the population aver- 



ages themselves, g(X), become random variables. To compute the expected (average) value 



of (151~1 [53]) . a 2 [n], one has to average over the distribution of g(X), which depends on the 
distribution /(A) of the A^s (we are assuming A-s i.i.d. random variables). This averages 
were obtained numerically, by evaluating expressions f|5T| |5U|) over the same realizations 
of the Aj's that were used in the numerical simulations. One can use the approximation 



g(X) ~ g(A), that works better the larger the N and the lower the variance a 2 , and that, 
due to the law of large numbers, is valid in the limit N — > oo. In Fig. 2 of the main text we 
compare the average of the analytical expression (I53I) with results coming from numerical 
simulations. We find perfect agreement and see that at first order the dependence of a 2 [n] 
with a 2 = A 2 — A 2 is linear and independent of the form of the distribution, as indicated by 



( I5T1) . Higher order corrections are noticeable for higher levels of diversity. 

In the case of heterogeneity in the preference of the agents for the states, as indicated in 
the main text, the variance is given by: 
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Fa + £ 
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a 2 \e] 



2e + A 



(54) 



In this case, the average of (15411 over the distribution of parameters can be easily computed, 
giving: 



a 2 [n] 
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St 



2 + X - 

e 



X_ 

2T 2 



2e + X/N X 



e(2e + A) 2e 2 A^ 



(55) 



The correlation function can be derived as follows (we exemplify the derivation in the 
case of distributed influence, for other types of heterogeneity, the derivation is similar): 
(145]) is an equation for the conditional averages (sj|{sz(to)}) if we set {si(t )} as initial 
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conditions. It implies: 



= eA - 2ea -> a(t + 1) = -(1 - e~ 2et ) + a(t )e~ 2et , 
at 2 



(56) 



with a = $2 fc Afc(sjfc|{si(t )})/^ r . Noticing that gSD is equal to ^ = e - (2e + X)(s l ) + a(t), 
we obtain: 



<*(*„ + t)\{s k (t )}} = \{l - e~^) + a(t0) - V2 e- 2 "(1 



e- A ') + ^(t )e- (2e+A)t . (57) 



Using nowK s( [n](f) = ((n(t +t)\n(t Q ))n(t )) st - (n) 2 t = E^((^(*o+OI{^(*o)})si(*o)) - x 
(remember (n)^ = N/2), and after some straightforward algebra, we obtain: 

K st [n](t) = (a 2 st - C/\)e-^ +J » + C/Xe~ 2 ^. (58) 



From (jUJ) we get C/X 



2e+X 



2XQ.-1/N) 

expression displayed in the main text. 



<y 2 t — N/4) = u, showing that (1381 is equal to the 



E. Appendix 



We start with equation ((4. 
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<k N 



(59) 

2(2e + A) V ; 

Using the rescaled variables = 4(7jj,Afc = . ttw ; an d defining 5^ := X]f?=o Aj Xj&ij, 



we obtain: 



2(2e+A)N' 

NX - 1 _ AT ^—7 

S'n+i = S n + y (A"Si + A"+! 

Defining now G n := ( -J— ) 5 n , T M := £) n=1 G n , we arrive to: 
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(61) 
(62) 



If limAf^oo G A / = 0, we see that: 
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TV v^oo / 2 
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(63) 
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Going back to the original variables, we finally obtain, with the notation of the main text: 

d= — - — , {,j+ rJ — . (64) 

2 Z^n=l \(\ + 4e)N J 



which can be rewritten in the form (I52j) . completing the proof 
The condition of convergence is: 

m'Too m'Too 2^ \ (A + 4e)iV / (2e + \)N 



lim G M = lim V _ 1 -^^cn j = 0. (65) 

■*-Hx>^\(\ + 4e)Nj (2e + A) AT 1 ; 



A necessary and sufficient condition for this is A, < ( A +^ e ) JV ; \/i — l ; ...iV. When the 
parameters A, are i.i.d. r. v. the probability of this typically approaches 1 as N grows. 
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